panel.grid.minor = element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x = element_text(angle = 70,hjust=.9,family="Times New Roman", color="black"),
axis.text.y = element_text(family="Times New Roman"),
axis.title.x = element_text(family="Times New Roman", size = 12),
axis.title.y = element_text(family="Times New Roman", size = 12))+
annotate(geom="text", x=1.5, y=98, label="Chi-square",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman")+ #adds text on xy coordinates
annotate(geom="text", x=1.5, y=92, label="p = 0.3508",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman") #adds text on xy coordinates
#stack the data to plot it
time=c("no travel ban","no travel ban","no travel ban","no travel ban","no travel ban","travelban","travelban","travelban","travelban","travelban")
Region= c("al-Hasakeh","Dara'a","Deir Ez-Zor","Hama","Homs","al-Hasakeh","Dara'a","Deir Ez-Zor","Hama","Homs")
Events=c(11,1,80,3,16,5,3,51,2,5)
tbanplot=data.frame(Region,time,Events)
tbanplot=tbanplot%>%
arrange(Region)
ggplot(data=tbanplot, aes(x=Region, y=Events, fill=time)) +
geom_bar(stat="identity", position=position_dodge())+
theme_bw()+
scale_fill_manual(values=c("#ace5ee", "red"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x = element_text(angle = 70,hjust=.9,family="Times New Roman", color="black"),
axis.text.y = element_text(family="Times New Roman"),
axis.title.x = element_text(family="Times New Roman", size = 12),
axis.title.y = element_text(family="Times New Roman", size = 12))+
annotate(geom="text", x=1.5, y=98, label="Chi-square",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman")+ #adds text on xy coordinates
annotate(geom="text", x=1.5, y=92, label="p = 0.3508",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman") #adds text on xy coordinates
rm(list=ls()) #Clearing workspace
#CURFEW
#Syria comparing numbers for before travelban and curfew (weeks 13-26 ) and during (weeks 65-78)
#data come from file supplied by Dawn 1.4.20
x=matrix(c(11,4,15,1,80,3,16,2,5,0,4,4,51,2,6,0), nrow=8)
chisq.test(x) #Chi-squared approximation may be incorrect
chisq.test(x, simulate.p.value = TRUE) #p-value = 0.4743
#stack the data to plot it
time=c("no curfew","no curfew","no curfew","no curfew","no curfew","no curfew","no curfew","no curfew","curfew","curfew","curfew","curfew","curfew","curfew","curfew","curfew")
Region= c("al-Hasakeh","Aleppo","Ar-Raqqa","Dara'a","Deir Ez-Zor","Hama","Homs","Idlib","al-Hasakeh","Aleppo","Ar-Raqqa","Dara'a","Deir Ez-Zor","Hama","Homs","Idlib")
Events=c(11,4,15,1,80,3,16,2,5,0,4,4,51,2,6,0)
tbanplot=data.frame(Region,time,Events)
tbanplot=tbanplot%>%
arrange(Region)
ggplot(data=tbanplot, aes(x=Region, y=Events, fill=time)) +
geom_bar(stat="identity", position=position_dodge())+
theme_bw()+
scale_fill_manual(values=c("#ace5ee", "red"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x = element_text(angle = 70,hjust=.9,family="Times New Roman", color="black"),
axis.text.y = element_text(family="Times New Roman"),
axis.title.x = element_text(family="Times New Roman", size = 12),
axis.title.y = element_text(family="Times New Roman", size = 12))+
annotate(geom="text", x=1.5, y=98, label="Chi-square",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman")+ #adds text on xy coordinates
annotate(geom="text", x=1.5, y=92, label="p = 0.1409",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman") #adds text on xy coordinates
rm(list=ls()) #Clearing workspace
#TRAVELBAN
#stack the data to plot it
time=c("no travelban","no travelban","travelban","travelban")
Region= c("North Sinai","Suez","North Sinai","Suez")
Events=c(69,1,74,0)
tbanplot=data.frame(Region,time,Events)
tbanplot=tbanplot%>%
arrange(Region)
ggplot(data=tbanplot, aes(x=Region, y=Events, fill=time)) +
geom_bar(stat="identity", position=position_dodge())+
theme_bw()+
scale_fill_manual(values=c("#ace5ee", "red"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x = element_text(angle = 70,hjust=.9,family="Times New Roman", color="black"),
axis.text.y = element_text(family="Times New Roman"),
axis.title.x = element_text(family="Times New Roman", size = 12),
axis.title.y = element_text(family="Times New Roman", size = 12))+
annotate(geom="text", x=1, y=98, label="Chi-square",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman")+ #adds text on xy coordinates
annotate(geom="text", x=1, y=92, label="p = 0.4743",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman") #adds text on xy coordinates
#CURFEW
#stack the data to plot it
time=c("no curfew","no curfew","curfew","curfew")
Region= c("North Sinai","Suez","North Sinai","Suez")
Events=c(69,1,74,0)
tbanplot=data.frame(Region,time,Events)
tbanplot=tbanplot%>%
arrange(Region)
ggplot(data=tbanplot, aes(x=Region, y=Events, fill=time)) +
geom_bar(stat="identity", position=position_dodge())+
theme_bw()+
scale_fill_manual(values=c("#ace5ee", "red"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x = element_text(angle = 70,hjust=.9,family="Times New Roman", color="black"),
axis.text.y = element_text(family="Times New Roman"),
axis.title.x = element_text(family="Times New Roman", size = 12),
axis.title.y = element_text(family="Times New Roman", size = 12))+
annotate(geom="text", x=1, y=98, label="Chi-square",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman")+ #adds text on xy coordinates
annotate(geom="text", x=1, y=92, label="p = 0.4743",
color="black", angle=0, size=4, fontface="bold",family="Times New Roman") #adds text on xy coordinates
install.packages("poliscidata")
#pg 15
#this line of code "switches on" the package so you can use it during this
#coding session. A library has to be installed once, but turned on everytime
#you start a new R session.
library(poliscidata)
#p15
#this does a fun little display welcoming you to the poliscidata packagee.
#feel free to skip it.
welcome()
# R Code
# Jóhanna K Birnir
# Objective: drawing functions from chapter 3
# Ernesto's Examples
############################################################################
#First example
c<-.5
p<-.3
op <- par(mfcol=c(1,2))
curve(p*log(-x^3+c/x+1), from=0, to=.33333, xlim=c(0,1), ylim=c(-1,1))
curve((p-1)*sqrt(abs(x-c/x)), from=.33, to=1, xlim=c(0,1), add=TRUE, col=2)
grid()
c<-.1
p<-.6
curve(p*log(-x^3+c/x+1), from=0, to=.33333, xlim=c(0,1), ylim=c(-1,1))
curve((p-1)*sqrt(abs(x-c/x)), from=.33, to=1, xlim=c(0,1), add=TRUE, col=2)
grid()
par(op)
c<-.5
p<-.3
curve(p*log(-x^3+c/x+1), from=0, to=.33333, xlim=c(0,1), ylim=c(-1,1))
curve((p-1)*sqrt(abs(x-c/x)), from=.33, to=1, xlim=c(0,1), add=TRUE, col=2)
c<-.1
p<-.6
curve(p*log(-x^3+c/x+1), from=0, to=.33333, xlim=c(0,1), ylim=c(-1,1), add=TRUE, lty=2)
curve((p-1)*sqrt(abs(x-c/x)), from=.33, to=1, xlim=c(0,1), add=TRUE, col=2, lty=2)
grid()
#CrossCutting function
p=.5
#Cross-Cutting function 1 with p
curve((p)*log10(-x^3+((1-x)/(2*x))+1), from=0, to=.3333, xlim=c(0,1), ylim=c(-.5,1), lty=1, xlab="Population balance: Overlapping cleavage", ylab="Minority utility")
#Cross-Cutting function 2 with 1-p
curve(-(1-p)*sqrt(abs(x^3-((1-x)/(2*x)))), from=.3333, to=1, xlim=c(0,1), ylim=c(-.5,1), col=1, lty=1, add=TRUE)
grid()
#these are the correct derivatives but they are not informative and so are omitted.
#Cross-Cutting derivative 1 with p not 1-p
curve(((p)*(2*3*(x^4)+1))/((2*x^5)-(x^2)-x), from=0, to=1, xlim=c(0,1), ylim=c(-1,1), col=2, lty=5, add=TRUE)
#Cross-Cutting derivative 2 with 1-p  CHECKED 12.17.18
curve(((1-p)*(2*(x^4)+1))/((4*(x^2))*((abs((x^3)-((1-x)/2*x)))^(1/2))), from=.0, to=1, xlim=c(0,1), ylim=c(-1,1), col=2, lty=5, add=TRUE)
#these are the correct derivatives but they are not informative and so are omitted.
#Cross-Cutting derivative 1 with p not 1-p
curve(((p)*(2*3*(x^4)+1))/((2*x^5)-(x^2)-x), from=0, to=1, xlim=c(0,1), ylim=c(-1,1), col=2, lty=5, add=TRUE)
#CrossCutting function
p=.5
#Cross-Cutting function 1 with p
curve((p)*log10(-x^3+((1-x)/(2*x))+1), from=0, to=.3333, xlim=c(0,1), ylim=c(-.5,1), lty=1, xlab="Population balance", ylab="Minority utility")
#Cross-Cutting function 2 with 1-p
curve(-(1-p)*sqrt(abs(x^3-((1-x)/(2*x)))), from=.3333, to=1, xlim=c(0,1), ylim=c(-.5,1), col=1, lty=1, add=TRUE)
#Cross-Cutting function 1 with p
curve((p)*log10(-x^3+((1-x)/(2*x))+1), from=0, to=.3333, xlim=c(0,1), ylim=c(-.5,1), lty=1, xlab="Population balance", ylab="Minority utility")
#Cross-Cutting function 2 with 1-p
curve(-(1-p)*sqrt(abs(x^3-((1-x)/(2*x)))), from=.3333, to=1, xlim=c(0,1), ylim=c(-.5,1), col=1, lty=1, add=TRUE)
grid()
#pg 15
#this line of code "switches on" the package so you can use it during this
#coding session. A library has to be installed once, but turned on everytime
#you start a new R session.
library(poliscidata)
freq(nes$abortlaw10)
#pg 20
#let's do another frequency table, this time with a graph
#the graph is generated automatically with the command
freq(nes$pid_x, nes$wt)
freq(nes$abortlaw10, nes$wt)
#2.16.21  THIS FILE IS NOT READY TO BE DISTRIBUTED.  I WILL REMOVE THIS COMMENT WHEN IT IS
# Code written by Johanna K. Birnir for Chapter 4
#This file replicates data composition and coding from earlier file
#### SETUP ####
library("tidyverse")
library("poliscidata")
library("readxl")
library("foreign")
options(tibble.print_max = Inf) #for rows
options(tibble.width = Inf) #for columns
rm(list=ls()) #Clearing workspace
#to be changed when releasing
maindir="~/Dropbox/Birnir_and_Satana/Alternatives_book/Cambridge/Chapters/Ch.4.Cross_national/Final.analysis.11.1.20/Overos.replication/cleaned.for.distribution.2.16.21"
setwd(paste0(maindir,"/original_data/"))  #where to get the data from
#svensson/relac
svensson18=read_excel("SvenssonandNilsson18/original__Svensson_18_data/Relac-JCRrep.xlsx")
getwd()
#svensson/relac
SaN18=read_excel("SvenssonandNilsson18/Relac-JCRrep.xlsx")
setwd(paste0(maindir,"/original_data/"))  #where to get the data from
#to be changed when releasing
maindir="~/Dropbox/Birnir_and_Satana/Alternatives_book/Cambridge/Chapters/Ch.4.Cross_national/Final.analysis.11.1.20/Overos.replication/cleaned.for.distribution.2.16.21"
setwd(paste0(maindir,"/original_data/"))  #where to get the data from
maindir="~/Dropbox/Birnir_and_Satana/Alternatives_book/Cambridge/Chapters/Ch.4.Cross_national/Final.analysis.11.1.20/Overos.replication/overos_ch4_rep11.4.20"
setwd(paste0(maindir,"/original data/"))  #where to get the data from
rm(list=ls()) #Clearing workspace
maindir="~/Dropbox/Birnir_and_Satana/Alternatives_book/Cambridge/Chapters/Ch.4.Cross_national/Final.analysis.11.1.20/Overos.replication/overos_ch4_rep11.4.20"
setwd(paste0(maindir,"/original data/"))  #where to get the data from
#Also possible to use latex tikz
library(DiagrammeR)
#name the objects and then specify which object has an arrow pointing to another object
grViz("digraph 'Data Merging Process' {
concentrate=true;
'relac \n (2812 org years)';
'amar \n (1202 groups)';
'relac \n (2812 org years)' -> 'relac.acd \n (2349 ethnic org years)';
'acd \n (518 orgs, 208 groups)' -> 'relac.acd \n (2349 ethnic org years)';
'acd \n (518 orgs, 208 groups)' -> 'amar.acd \n (217 groups)';
'amar \n (1202 groups)' -> 'amar.acd \n (217 groups)';
'amar.acd \n (217 groups)' -> 'relac.acd.amar \n (2016 group years)';
'relac.acd \n (2349 ethnic org years)' -> 'relac.acd.amar \n (2016 group years)';
'relac.acd.amar \n (2016 group years)' -> 'relac.acd.amar.areligion \n (1970 group years)';
'a-religion \n (1199 groups)' -> 'relac.acd.amar.areligion \n (1970 group years)';
'relac.acd.amar.areligion \n (1970 group years)' -> 'final.df \n (1970 group years)';
'Sambanis, Penn, Jones & Lupu, Selway \n (various)' -> controls;
controls -> 'final.df \n (1970 group years)';
labelloc='t';
label='The Data Merging Process: \n Number and Unit in Parentheses';
}
")
#name the objects and then specify which object has an arrow pointing to another object
grViz("digraph 'Data Merging Process' {
concentrate=true;
'RELAC-domestic \n (420 orgs, 2044 years)';
'ACD2 \n (321 orgs, 208 groups)';
'A-Religion \n (1199 groups)';
'RELAC-domestic \n (420 orgs, 2044 years)' -> 'RELAC.ACD2 \n (207 groups, 2016 years)';
'ACD2 \n (321 orgs, 208 groups)' -> 'RELAC.ACD2 \n (207 groups, 2016 years)';
'ACD2 \n (321 orgs, 208 groups)' -> 'A-Religion.ACD2 \n (220 groups)';
'A-Religion \n (1199 groups)' -> 'A-Religion.ACD2 \n (220 groups)';
'A-Religion.ACD2 \n (220 groups)' -> 'RELAC.A-Religion \n (1970 group years)';
'RELAC.ACD2 \n (207 groups, 2016 years)' -> 'RELAC.A-Religion \n (1970 group years)';
'RELAC.A-Religion \n (1970 group years)' -> 'final.df';
'Penn, Selway, Jones & Lupu,  Sambanis \n (various)' -> 'controls \n (1273 country years)';
'controls \n (1273 country years)' -> 'final.df';
{rank=same;'RELAC-domestic \n (420 orgs, 2044 years)'; 'ACD2 \n (321 orgs, 208 groups)'; 'A-Religion \n (1199 groups)'}
}
")
#name the objects and then specify which object has an arrow pointing to another object
grViz("digraph 'Data Merging Process' {
concentrate=true;
'RELAC-domestic \n (420 orgs, 2044 years)';
'A-Religion \n (1199 groups)';
'RELAC.A-Religion \n (201 groups, 1970 group years)' -> 'final.df';
'Penn, Selway, Jones & Lupu,  Sambanis \n (various)' -> 'controls \n (1273 country years)';
'controls \n (1273 country years)' -> 'final.df';
{rank=same;'RELAC-domestic \n (420 orgs, 2044 years)'; 'A-Religion \n (1199 groups)'}
}
")
#name the objects and then specify which object has an arrow pointing to another object
grViz("digraph 'Data Merging Process' {
concentrate=true;
'RELAC-domestic \n (420 orgs, 2044 years)' -> 'RELAC.A-Religion \n (201 groups, 1970 group years)';
'A-Religion \n (1199 groups)'-> 'RELAC.A-Religion \n (201 groups, 1970 group years)';
'RELAC.A-Religion \n (201 groups, 1970 group years)' -> 'final.df';
'Penn, Selway, Jones & Lupu,  Sambanis \n (various)' -> 'controls \n (1273 country years)';
'controls \n (1273 country years)' -> 'final.df';
{rank=same;'RELAC-domestic \n (420 orgs, 2044 years)'; 'A-Religion \n (1199 groups)'}
}
")
rm(list=objects()) #clear the environment
# Set your working directory
setwd("/Users/folder")
library(poliscidata)
wtd.t.test(gss$egalit_scale, 0, weight=gss$wtss)
CI95(19.3445353, 0.3039164)
wtd.t.test(gss$egalit_scale, 20.1, weight=gss$wtss)
CI95(-0.7554647, 0.3039164) #the 95% CI for the mean difference of egalitarian scores between typical adult and social science majors
wtd.t.test(gss$egalit_scale, 18.8, weight=gss$wtss)
CI95(0.5445353, 0.3039164)
#Q2. A
wtd.ttestC(~int_info_scale, ~age2, gssD)
wtd.t.test(gss$egalit_scale, 0, weight=gss$wtss)
CI95(19.3445353, 0.3039164)
wtd.t.test(gss$egalit_scale, 20.1, weight=gss$wtss)
install.package(haven)
install.packages("haven")
#7.27.21
# Code written by Johanna K. Birnir
#Purpose: read and explore CIP sav file
#### SETUP ####
library("tidyverse")
library("haven)" #to read .sac
install.packages("readspss")
install.packages(readspss)
install.packages('readspss')
library(poliscidata)                   # Loads R companion package in session
# CROSS-TABULATION ANALYSIS WITH A CONTROL VARIABLE
#pg 77
# Cross-tabulation without control variable
#legalize marijuana by religious service attendance
xtp (gss, grass, attend3, wtss)
# Cross-tabulation with control variable
#syntax xtabC(~DV + IV + CV, design.dataet)
#pg 78
xtabC (~grass + attend3 + kids, gssD)
#move across top row for kids no, drops about 34 points.
#move across top row for kids yes, drops about 31 points
#so, roughly the same
#pg 79
#Controlled cross-tabulation, interaction relationship
xtabC (~divlaw2 + attend3 + kids, gssD)
# MULTIPLE LINE CHARTS
#pg 80
# Retrieve level names of variable to graph.
levels(gss$grass)
#Create the indicator in the non-design dataset, multiplied times 100.
#this line pulls out the observations that have "LEGAL" as their value on this variable
#and gives them a value of 1
gss$grass.legal=as.numeric(gss$grass=="LEGAL")
#make it a percent
gss$grass.legal=100*gss$grass.legal
#Run svydesign, updating gssD to include the new indicator.
gssD=svydesign(id=~1,data=gss,weights=~wtss)
# Multiple line chart, controlled comparison
#syntax iplotC(~DV, ~IV+ CV, design.dataset, DV~CV+IV, plot options)
#NOTE HOW THE IV AND CV SWITCH PLACES IN THE COMMAND
iplotC(~grass.legal,~attend3+kids,gssD,grass.legal~kids+attend3,
xlab="Relgious attendance",
ylab="Percent favoring legalization",
main=list("Percentage Favoring Marijuana Legalization\n by Religious
Attendance and Kids",cex=1))
# THE legend FUNCTION
#pg 81
# the code below layers on top of line chart created above
legend("topright",   #Locates the legend in the plotting area
legend=c("No","Yes"),#Labels the categories
lty=c(1,2), 	   #Make sure the lin types match the legend labels
lwd=2, 		   #Specifies the same line weight shown in the graph
title="Does R have kids?",
inset=0.1, 	   #Specifies distance from the plot borders
bty="n")          #Suppresses the box around the legend
# MEAN COMPARISON ANALYSIS WITH A CONTROL VARIABLE
#p 82
# Controlled mean comparison
#syntax imeansC(~DV, ~IV+CV, design.dataset)
imeansC (~ft_dem, ~gender + married, nesD)
#pg 83
#do the plot
iplotC(~ft_dem,~gender+married,nesD,ft_dem~married+gender,
xlab="Gender", ylab="Democratic Party Rating",
main="Democratic Party Ratings\nby Gender and Marital Status")
legend("topleft",legend=c("No","Yes"),lty=c(1,2),lwd=2,
title="Is R married?",inset=0.1, bty="n")
#pg 84
#get the means
imeansC(~modsex_scale, ~link_wom_scale + dem_raceeth2, nesD)
#pg 85
#do the plot
iplotC(~modsex_scale,~link_wom_scale+dem_raceeth2,nesD,
modsex_scale~dem_raceeth2+link_wom_scale,
xlab="Linked Fate Scale", ylab="Modern Sexism Scale",
main=list("Modern Sexism Opinions by \n Linked Fate and Race",cex=1.2))
legend("topleft",legend=c("Black","White"),lty=c(2,1),lwd=2,
title="Race",inset=0.1, bty="n")
load("~/Downloads/pilot.data (1).Rdata")
names(pilot.data)
knitr::opts_chunk$set(echo = TRUE)
summary(cars)
knitr::opts_chunk$set(echo = TRUE)
summary(cars)
plot(pressure)
plot(pressure)
plot(pressure)
plot(pressure)
summary(cars)
summary(cars)
plot(pressure)
summary(cars)
summary(cars)
knitr::opts_chunk$set(echo = TRUE)
summary(cars)
plot(pressure)
plot(pressure)
plot(pressure)
summary(cars)
knitr::opts_chunk$set(echo = TRUE)
summary(cars)
plot(pressure)
install.packages("blogdown")
library(blogdown)
load("/Users/hannahbirnir/Dropbox/Birnir_and_Satana/CWC_paper/CWC_JOP/df.joprnr12722.Rdata")
names(df.joprnr12722.Rdata)
names(df.joprnr12722)
df.joprnr12722=df.joprnr12722%>%
mutate(Uganda=ifelse(amar_country=="Uganda",1,0))
library("tidyverse")
library("stargazer")
library("sandwich")
library("lmtest")
library("margins")
df.joprnr12722=df.joprnr12722%>%
mutate(Uganda=ifelse(amar_country=="Uganda",1,0))
names(df.joprnr12722)
CWC_supplementary=df.joprnr12722%>%
select("year","Religious.incompatibility", "Shared.family.ns.nd",  "Population.balance.nod",
"Organizational.competition","National.cross.cutting",  "ln.GDP.per.capita",
"Democracy", "Separatism", "T1", "T2", "T3", "Shared.sect.nd", "Nepal", "Uganda")
load("/Users/hannahbirnir/Dropbox/Birnir_and_Satana/CWC_paper/CWC_JOP/replication.materials/CWC_replication/df.CWC.622.Rdata")
#Checking the data
model3 = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3,
family = "binomial", data = df.CWC.622)
model3.robust.SE=coeftest(model3, vcov = vcovHC(model3, sandwich=TRUE))
stargazer(model3.robust.SE, title="Results", type = "text",align=TRUE)
model3.check = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3,
family = "binomial", data = CWC_supplementary)
model3.robust.SE.check=coeftest(model3.check, vcov = vcovHC(model3.check, sandwich=TRUE))
stargazer(model3.robust.SE.check, title="Results", type = "text",align=TRUE)
model1s = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3,
family = "binomial", data = CWC_supplementary)
model1s.robust.SE=coeftest(model1s, vcov = vcovHC(model1s, sandwich=TRUE))
stargazer(model3.robust.SE.check, title="Results", type = "text",align=TRUE)
model1s = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3, Nepal, Uganda,
family = "binomial", data = CWC_supplementary)
model1s = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3+Nepal+Uganda,
family = "binomial", data = CWC_supplementary)
model1s.robust.SE=coeftest(model1s, vcov = vcovHC(model1s, sandwich=TRUE))
stargazer(model3.robust.SE.check, title="Results", type = "text",align=TRUE)
stargazer(model1s.robust.SE.check, title="Results", type = "text",align=TRUE)
model1s = glm(Religious.incompatibility ~ Shared.family.ns.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3+Nepal+Uganda,
family = "binomial", data = CWC_supplementary)
model1s.robust.SE=coeftest(model1s, vcov = vcovHC(model1s, sandwich=TRUE))
stargazer(model1s.robust.SE.check, title="Results", type = "text",align=TRUE)
stargazer(model1s.robust.SE, title="Results", type = "text",align=TRUE)
model2s = glm(Religious.incompatibility ~ Shared.sect.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3,
family = "binomial", data = CWC_supplementary)
model2s.robust.SE=coeftest(model1s, vcov = vcovHC(model2s, sandwich=TRUE))
stargazer(model2s.robust.SE, title="Results", type = "text",align=TRUE)
names(CWC_supplementary)
model2s = glm(Religious.incompatibility ~ Shared.sect.nd*Population.balance.nod+Organizational.competition+
National.cross.cutting+ln.GDP.per.capita+Democracy+Separatism+T1+T2+T3,
family = "binomial", data = CWC_supplementary)
model2s.robust.SE=coeftest(model2s, vcov = vcovHC(model2s, sandwich=TRUE))
stargazer(model2s.robust.SE, title="Results", type = "text",align=TRUE)
load("/Users/hannahbirnir/Dropbox/Birnir_and_Satana/CWC_paper/CWC_JOP/df.joprnr12722.Rdata")
CWC_supplementary_cross=df.joprnr12722%>%
select("year","Religious.incompatibility", "Shared.family.ns.nd",  "Population.balance.nod", "numcode")%>%
arrange(numcode, year)
CWC_supplementary_cross=CWC_supplementary_cross%>%
group_by(numcode)%>%
mutate(id = seq_along(numcode))
CWC_supplementary_cross=CWC_supplementary_cross%>%
group_by(numcode)%>%
mutate(id = seq_along(numcode))
View(CWC_supplementary_cross)
CWC_supplementary_cross=CWC_supplementary_cross%>%
group_by(numcode)%>%
mutate(id = seq_along(numcode))%>%
filter(id==1)
#Locking down violence 10.18.22
#Replication file
#Appendix join-count
####################################################################################
#Packages
library(tidyverse)
library("spdep")
library("raster")
load(file= "nolockdown.spdf", verbose = FALSE)
nolockdown.spdf <- nolockdown
coordinates(nolockdown.spdf) <- c("longitude", "latitude")
load(file= "lockdown.spdf", verbose = FALSE)
lockdown.spdf <- lockdown
coordinates(lockdown.spdf) <- c("longitude", "latitude")
setwd("~/Dropbox/COVID19peace/APSR/Replication Files (APSR)/GIS/GIS_appendix")
load(file= "nolockdown.spdf", verbose = FALSE)
nolockdown.spdf <- nolockdown
coordinates(nolockdown.spdf) <- c("longitude", "latitude")
load(file= "lockdown.spdf", verbose = FALSE)
lockdown.spdf <- lockdown
coordinates(lockdown.spdf) <- c("longitude", "latitude")
#Define the extent for the join count analyses for lockdown
jc.extent <- extent(39,46,30,40)
#Set cell size
r <- raster(nrows=25, ncols=20, ext=jc.extent)
#Join Count Analysis for NO LOCKDOWN
nolockdown.rast <- rasterize(nolockdown.spdf, r, field = 1)
nolockdown.rast[is.na(nolockdown.rast)] <- 0
nbnl <- cell2nb(nrow = nrow(nolockdown.rast), ncol = ncol(nolockdown.rast))
# Convert the neighbors list to a 'weights' list;
lwbnl <- nb2listw(nbnl, style = "B")
joincount.test(as.factor(nolockdown.rast@data@values), lwbnl, alternative = "greater")
#Join Count Analysis for LOCKDOWN
lockdown.rast <- rasterize(lockdown.spdf, r, field = 1)
lockdown.rast[is.na(lockdown.rast)] <- 0
nbnl <- cell2nb(nrow = nrow(lockdown.rast), ncol = ncol(lockdown.rast))
# Convert the neighbors list to a 'weights' list;
lwbnl <- nb2listw(nbnl, style = "B")
joincount.test(as.factor(lockdown.rast@data@values), lwbnl, alternative = "greater")
